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ABSTRACT 



Landing aircraft on board carriers is a most delicate phase of flight operations at 
sea. The ability to predict the aircraft carrier's motion over an interval of several sec- 
onds within reasonable error bounds may allow improvement in touchdown dispersion 
and a more certain value for a ramp clearance due to a smoother aircraft trajectory. 
Also, improved information to the Landing Signal Officer should decrease the number 
of waveoffs. 

This work indicates and shows graphically that, based on the data for pitch, heave 
and roll measured for various ships and sea conditons, the motion can be predicted well. 
The predictor was designed on the basis of Kalman's optimum filtering theory for the 
discrete time case, adapted for real-time digital computer operation. 
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I. INTRODUCTION 



The landing phase of VTOL, STOL and HELO aircraft aboard a ship without a 
ship's Inertial Navigational System represents a complex operation and a demanding 
task. The last 10 to 15 seconds before the aircraft touchdown involves terminal guidance 
and control problems, where not only the aircraft is disturbed by several kinds of 
stochastic disturbances (wind), but also the touchdown point (on the ship) is being 
moved randomly. Despite the wind disturbances and the final point (target) random 
motion, the landing accuracy specified for ship operations is very high, i.e., a few tens 
of feet longtitudinal landing dispersion. Such a terminal point problem is made tractable 
in a most natural way by assuming that the ship's position can be predicted for several 
seconds ahead so that the aircraft is guided toward the future position of the touchdown 
point. The scope of this study was to establish to what extent a stochastic process, like 
the ship's motion, is predictable over moderate periods of time. 

Graphical results obtained throughout this estimator's feasibility study concerning 
the relationship between the ship's states: heave, pitch and roll versus the ship's esti- 
mated states; heave, pitch and roll; and the influence of noise measurement are pre- 
sented. Digital simulations show that the prediction accuracy was very accurate for 
different swell wave height and period. The feasibility of estimating the ship's motion 
of acceptable bounds of error can also lead to improvement of the Landing Signal Offi- 
cer's (LSO) decision policy for waveoffs. 

This study was based on the use of the ship's motion equations and the simulated 
measuring instrumentations existing on board ship in order to get the predicted motions. 
Using this information, a predictor based on Kalman's theory of optimum estimation 
was designed. 

Several circumstances contributed to the success of this approach. The size and 
mass of the ship significantly filter the motion of the sea. A complete landing operation 
is short enough that the stochastic processes are reasonably taken to be stationary. 
Finally, the prediction interval is only a small fraction of the time it takes each aircraft 
for final approach. 

This work is divided into three parts: 1) the derivation of the mathematical model 
of the ship's motion, 2) the rationale in the implementation of Kalman filter and pre- 
dictor equations, and 3) discussion of the simulated computer results obtained. Since 
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we are interested only in the most critical aspects of the landing operation, namely the 
characteristics of the longtitudinal channel, we merely investigate in the sequel the pre- 
dictability of the heave, pitch and roll motion of the ship where yaw is omitted. 
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II. PROBLEM STATEMENT 



A. GENERAL 

The ship's attitude estimation used in this thesis involves calculation of heave, pitch 
and roll of the ship moving in the general direction towards a swell while recovering an 
aircraft. The motion of the ship is given in z ,6 and 0 coordinates as shown in Figure 

This problem will be developed using state space methods. Given the heave, pitch 
and roll motion (the measurements) received through the ship's sensor system, we are 
interested in estimating the heave, pitch and roll of the ship. The state variables for this 
plant are z h , i A , 6 f , 8 f , 0, and 0, . 

B. SHIP MOTION 

A ship moving on the surface of the sea is almost always in oscillatory motion. The 
different kinds of oscillatory motions that a ship experiences can be described with the 
help of Figure 1, which shows the six kinds of motion, three linear and three rotational 
about the principal axes. Accordingly, 

The six kinds of ship's motions are: 

1. a = surging - motion backwards in the direction of the ship travel. 

2. b = staying -* athwartship motion of the ship. 

3. c = heaving — motion vertically up and down. 

4. d = rolling -- angular motion about the longitudinal axis. When the ship rolls, it 
lists alternately from starboard to port and then back to starboard. 

5. e = pitching — angular motion about the transverse axis. When a ship pitches, it 
trims alternately by the bow and by the stern. 

6. / = yawing -- angular motion about the vertical axis. 
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x-axis is longitudinal axis 
j'-axis is transverse axis 
z-axis is vertical axis 



Figure 1. The x, y, and z-axes of a ship 
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Only three kinds of motion, namely heaving, pitching and rolling ( z , 6 , <£ ) are 
purely oscillatory motions, since these motions act under a restoring force or moment 
when the ship is disturbed from its equilibrium position. These ship motions are the only 
motions being considered in this thesis. 

Although in reality a ship experiences all six kinds of motion simultaneously, only 
one motion of the three angular motions will be treated at a time in the following 
sections. However, one must bear in mind that any one kind of motion is not inde- 
pendent of the others; in the first approximation, coupling between the motions is neg- 
lected in order to simplify the problem. 

C. SYSTEM MODEL 

The system to be modeled in this problem is the motion of a surface ship at sea 
during an aircraft recovery phase. This is a linear, time-angle system that can be de- 
scribed in a discrete equation of ship motion . The state space equation is 

&k + 1 = "b k u k '(2.1) 



where 

x k = parameter to be estimated (state vector) heave, pitch and roll 

d>*= state transition matrix which describes how the states of the dynamic system are 
related. 

T*= system random uncorrelated coefficient matrix. 
u k = sequence of random uncorrelated inputs. 

From Equation (2.1) and the above assumptions, the state vector for heave, pitch and 
roll is 



*k = 



6 

0 



4 > 






(2.2) 



and the system state equation can be approximated for T < < 1 
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where 

u k = £ a cos co s t + plant noise 

C„= swell wave amplitude 

co = swell angular frequency 

c, = g 3 = heaving restoring force coefficient 

b 2 = heaving damping force coefficient 

c p = gp~ pitching restoring moment coefficient 

b„ — pitching damping moment coefficient 

c, = g= rolling restoring moment coefficient 

b,= heaving damping moment coefficient. 

The system noise process for the ship motion estimation problem is a function of the 
random uncorrelated coefficient matrix. T* , and the random forcing function, u k . 

D. MEASUREMENT MODEL 

For a linear measurement process, the measurements are linearly related to the state 
variables and can be modeled using the linear measurement equation 

ik = H k*k + Jfc (2-4) 



where 

z k - set of measurements. 

II k = observation matrix that gives the noiseless relationship between the measure- 
ments and the state vector. 

x k = state vector 

v k = sequence of random uncorrelated measurement noise. 

In this ship motion estimation problem, the measurements are the heave, pitch and 
roll angles taken from the ship sensors. 
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z/< = m oD^ + vj 



(2.5) 



The linear equation with the measurement data available can now be processed with 
the Kalman filter. The type of noise that afTect the ship motion data are assumed to be 
white gaussian noise. This noise model was used in the computer simulations and shown 
in Figure 2, with zero mean and variance of one. 

The complete development of this model can be found in Bhattacharyya [Ref. 1: 
pp.35-180], Wall [Ref. 2: pp.31-47], Korvin- Kroukovsky [Ref. 3] and Blagoveshchensky 
[Ref. 4). 

The corresponding block diagram of the system is shown in Figure 3. 



Random Noise vs Time 




Time (sec)xlO 



Figure 2. White Noise Model 
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Figure 3. Ship Motion Model Block Diagram 



III. KALMAN FILTER THEORY 



A. GENERAL 

Filtering refers to the process of estimating the state vector at the current time based 
upon all past measurements. An optimal filter concentrates on optimizing a specific 
performance measure used to determine the quality of the estimate. The Kalman filter 
is an optimal filter in a class of linear filters that minimize the mean square estimation 
error between the actual and desired output. In other words, the Kalman filter attempts 
to minimize the elements along the main diagonal of the state error covariance matrix. 
The filter is a recursive algorithm for processing discrete measurements or observations 
in an optimal manner, [Ref. 5: p. 101] and [Ref. 6]. It requires a priori knowledge of the 
state estimate (v*^) and its error covariance (P k _ ]), and the current observation ( z k ). The 
Kalman filter is the proper algorithm to be used when both the system model and the 
measurement model are linear functions of the state variables and these models can be 
described by the equations 



In this equation, and T* are constant matrices and H A is a linear constant func- 
tion of the state variable x k , while u k is the random forcing function and v k is the meas- 
urement noise. The plant random forcing function and measurement noise are assumed 
to be uncorrelated (white Gaussian) with zero mean. That is 



&k+\ = ®k*k + 



(3.1) 



Zk = H k*k + Xk 



(3.2) 



= Q\k)6 u 



(3.3) 



and 



£M*).v r (/)] = R'Whj 



(3.4) 



where 



hj= 1 , k-j 

= 0, k*J 



(3.5) 



9 



B. NOISE PROCESSES 

The calculation of the error covariance matrix and the filter gain matrix requires the 
covariance matrices for the uncorrelated noise process u k and v k . For the measurement 
noise process, v k , the covariance matrix is 



where R k is defined as the state measurement noise covariance matrix. It is based on the 
sensor accuracy and accounts for unknown disturbances such as steps, white noise, or 
imperfections in the plant model. The variance of the white noise model used in the 
computer simulations was one. 

The state excitation matrix, Q k , used in the Kalman filter represents the system noise 
process and is a function of the system noise coefficient matrix, T*, and the random 
forcing function, u k . This matrix is given by 



where T A is the same as in Equation (2.3). The Q k matrix allows for any random heave, 
pitch and roll ship motion well as inaccuracies in the system model. The magnitude of 
O k has a direct bearing on the magnitude of the state error covariance matrix and it 
prevents the covariance matrix from becoming singular by ensuring some uncertainty in 
the state estimates. 

C. INITIALIZATION AND OPERATION 

In the ship motion estimation, Kalman filters are used to minimize the ship's atti- 
tude errors. Prior to processing the measurement data, the filter must be initialized with 
an initial state estimate and an initial error covariance matrix. This initialization process 
is a very important step in the filter operation and gross inaccuracies in this step may 
cause the filter to diverge. Divergence occurs when the calculated covariance errors be- 
come much smaller than the actual covariance errors. This causes the estimate values 
of the states to pull away from the actual value and is also used as a startup value for 
the recursive scheme. 

The basic operation of the filter is a relatively straightforward recursive process. 
Based on Kirk [Ref. 6 : pp.4-1, 4-62] and Gelb [Ref. 5], the discrete time version of 
equations used in the Ka!:...m filter are: 



El v k v fl = Rk 



(3.6) 



Qk = U~kQ'iJ~kl 



(3.7) 



£(k\k—\) - ^k^ikW 



(3.8) 
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P(k\k-l) = ^kP{k\k)^k + Qk 

G k = P { k\k-i)H T k {H k P {k \ k -,)H T k + R k ) 

^ik\k) ~ ^</c|/c — t) G k(-k ~~ H]JC(k\k-l)) 

P{k\k ) = (I~ G kHk)P{k\k-\) 



(3.12) 



(3.10) 



(3.11) 



(3.9) 



where 

= projected ahead state estimate 
<t> k = state transition matrix given by Eq. (2.3) 

= projected ahead state error covariance matrix 
Q k = state excitation covariance matrix given by Eq. (3.4) 

G k = Kalman gain matrix 

R k = state measurement noise covariance matrix given by Eq. (3.3) 

H k = linearized measurement matrix given by Eq. (2.5). 

From the given equations, we can start the filter processing operation. The a priori 
state estimate are calculated using the O matrix shown in Equation (2.3) where T is the 
sampling time in seconds between the observed ship motions (heave, pitch and roll). 
(It is assumed ship motions are received simultaneously through the ship sensors.) 

The Kalman gain matrix serves to minimize the mean square estimation error and 
is an indication of how much emphasis or weight will be placed on the current observa- 
tion. If P {kk -D is small, the Kalman gain matrix will also be small due to the finite value 
of R k . If the is relatively large, the gain is approximately one. By rewriting the 

equation for the calculation of the state estimate Equation (3.11), as 



we can see how the Kalman gain matrix directly affects the weight placed on the current 
observation z k . A large gain, indicating a large a priori error covariance, will place more 
weight on the current observation as the filter tries to correct the states. A small gain, 
indicating a small error covariance, places less emphasis on the new observation. 

If the Kalman gain is expressed as 



£(k\k) ~ (1 - G kHk)^(k\k-\) + G k±k 



(3.13) 




(3.14) 
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it can be seen that the gain matrix is "proportional" to the uncertainty in the estimate 
P(k k -\) and "inversely proportional" to the measurement noise R k . For a large R k and a 
small P {kk . X) , the measurement in Equation (3.2) is due mainly to noise and only small 
corrections should be made in the state estimate. However, if R k is small and P m -x> is 
large, the measurement contains considerable information about the errors in the esti- 
mates and therefore, a strong correction should be made to the state estimates [Ref. 5: 
pp. 127-8]. 

The computation process for prediction is divided into the following steps: 1) cal- 
culate the estimate using Equation (3.1 1) and 2) use Equation (3.8) for the desired pre- 
diction. 

The block diagram of the system with Kalman Filter is shown in Figure 4 and the 
system predictor block diagram is shown in Figure 5. 
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figure 4. Ship Motion Model Diagram With Estimator 
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Figure 5. Ship's Motion Predictor Block Diagram 



14 



IV. COMPUTER SIMULATIONS 



A. GENERAL 

A measured simulated ship motion (heave, pitch and roll) and the forcing function 
data were chosen based from the ship speed of 25 knots, heading 000 degrees true and 
the swell w r ave heading direction at 120 degrees true. The sea state for a sw'ell w r ave 
height of 11 feet is between 5 and 6 based on Bhattacharyya [Ref. 1: p. 104] which can 
be described that large w-aves begins to form while white crests are extensive everywhere 
(probably some spray). The wind force (beaufort) is classified as level 6 w T hich is char- 
acterized by a strong breeze with the velocity of 25 knots. 

These measured data w’ere used in the Matlab computer program ([Ref. 7] and [Ref. 
8]) algorithm (as shown in Appendix A) to generate the graphical results to see the 
physical representation of the ship motion oscillations. These graphical representations 
of swell, heave, pitch and roll, wilich are the results from the above computer runs, have 
the following parameters listed in Table 1 and are shown in Figures 6, 7, 8 and 9. 



Table 1. WAVE AND SHIP MOTION PARAMETERS, TS = 0.1 SEC 



Coefficient 


Swell 


Heave 


Pitch 


Roll 


Wave ht (ft) 
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- 


Period (sec) 


7 


11 


10 


18 


c o (rad sec) 


0.8976 


0.5712 


0.62S3 


0.3491 


Damping 


0 


0.6854 


0.7540 


0.4189 


Restoring 


0 


0.5712 


0.6283 


0.3491 


Kalman Gain 


-- 


0.1255 

0.0841 


0.1191 

0.0756 


0.2197 

0.2721 
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With the simulated observed ship motion data for heave, pitch and roll generated 
from the above computer runs, the Kalman filter was used to calculate the estimates by 
using the Matlab computer program algorithm generated as shown in Appendix B. The 
resulting estimated data were represented in graphical form to show the difference be- 
tween the observed and the estimated ship motions. These graphs, with estimates, are 
shown in Figures 10, 11 and 12. Computer runs w r ere also conducted to show the 
graphical representations of the predicted ship motions heave, pitch and roll. A Matlab 
computer program is shown in Appendix C for this calculation. These graphical pred- 
ictions of ship motions are shown in Figures 13, 14 and 15. As we can see, the estimates 
and predictions illustrate a well-defined representation of the ship motion that can be 
used in the prediction of the ship motion in the future time which is to be relayed to the 
Landing Signal Officer for his decision for waveoffs and to the aircraft's pilot through 
the aircraft's TACAN during the recovery phase of flight operations. 
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Figure 12 



Observed and Estimated Roll Motion caused by 11 feet high of Swell 
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Ships Observed Pitch and Predicted Pitch 




Figure 14. Observed and Predicted Pilch Motion caused by 11 feet high of Swell 
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V. CONCLUSIONS 



The feasibility of estimating an aircraft carrier motion at sea by measuring the ship's 
actual position and motion was investigated. The ship's motion representation was 
generated based on ship's motion mathematical model taken from [Ref. 1]. Subse- 
quently, a Kalman filter-estimator adapted for real-time computation on a digital com- 
puter, was generated. The results obtained show that the desired prediction time can 
be reached with reasonably acceptable errors. 

Being able to predict the ship's motion can lead to an improvement of aircraft 
landing accuracy and safety. This can be accomplished, for instance, by generating ter- 
minal guidance (landing) laws making use of the future ship's motion and position. 
Moreover, the possibility of prediction can improve the LSO information and policy for 
landing acceptance or wave-offs. The possibility of processing the measured motion by 
Fast Fourier Transform Algorithms (FFT), in order to obtain the estimated ship motion 
parameters such as the heave ( z ), pitch (6) and roll (<f>), in real time may lead toward 
an adaptive predictor scheme. 
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APPENDIX A. MATLAB PROGRAM 



This Matlab computer program calculates the Observed Ship's Motion such as the 
heave, pitch and roll of the ship. Required coefficients must be available to run the 
program. 



% A Matlab program to generate the ship's heave, pitch and roll 
% roll motion with a given Swell wave force of N sin ws*t. 

% 

% 

%%%%%%%%%%%%%%%%%%%%%%%%%% SWELL parameter %%%%%%%%%%%%%%%%%%%%%%%%%%% 

% 

% 

% 

diary p. m 

% Swell period (Ts) seconds 
Ts=7; 



% Swell angular frequency 

ws= 2*pi*l/Ts ; 



7o Wave Height (N) feet 
N=ll; 
bl=N*ws; 

% Restoring coefficient 



%%%%%%%%%%%%%%%%%%%%%%%%%% Heave parameter %%%%%%%%%%%%%%%%%%%%%%%%%%% 

% 

% Heave period (th) seconds typical for carriers 
th= 11; 

% 

% Heave angular frequency (wh) rad/sec 
wh= 2*pi*l/th ; 

% 

% Heave restoring cefficient (cz) 
cz=wh 

% Heave damping factor (zetah) 
zetah= 0. 6; 

% Damping coefficient (bh) 

bz=2*zetah*wh ; 

% 

% 

%%%%%%%%%%%%%%%%%%%%%%%%% Pitch parameter %%%%%%%%%%%%%%%%% %%%%%%%%%% 

% 
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% 

% Pitch period (tp) seconds typical for carriers 
tp=10; 



% Pitch angular frequency (wp) rad/sec 

wp= 2*pi*l/tp ; 

% 

% 

% Pitch Restoring coefficient 

cp=wp 

% J 

% Heave damping factor (zetap) 

zetap= 0. 6 ; 

% 

% Damping coefficient (bp) 

bp=2*zetap*wp ; 

% 

% 

%%%%%%%%%%%%%%%%%%%%%%%%% Rol 1 parameter %%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% 

% 

% Roll period (tr) seconds typical for carriers 
tr= 18; 

7 a Roll angular frequency (wr) rad/sec 
wr= 2*pi*l/tr ; 



% Roll Restoring coefficient 
cr=wr 
7a 

% Roll damping factor (zetar) 
zetar= 0. 6 ; 

7a Damping coefficient (br) 
br=2*zetar*wr 



7a 



7aVL7a7a7a7aV/a7a7aV/a7a7a7a7a7a7a7a7a7a7a Parameters Of the system 7a7a7aV/a7a7aV/a7a7a7a7a7aV/a7aT/a7a 
7a 

7a where: 

% a= []‘> -- the A martix coefficient in continous time 

% b=[];-- the B matrix coefficient in continuos time 

7a 

7a Swell Parameters 

7a ======== 



7a 

% 

7a 

7a 

7a 



as=-0 1; -cs 0-; 
bs=[ 0; 1] ; 



Heave Parameters 
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X 

% 

X 

% 



X 

X 

X 

% 

% 



diary 

% 

X 

% 

% 

% 



ah=[ 0 1; -02**2 -bz] ; 
bh=[ 0; . 3] ; 

Pitch Parameters 



ap=[0 1; -cp**2 -bp] ; 
bp=[ 0; . 3] ; 



Roll Parameters 



ar=[ 0 l;-cr**2 -br] ; 
br=[ 0; . 6] ; 
off 



Measurement Parameters 



:=[ 1 0 ] 
HO] 



dt=0. 1; 
q=eye( 2); 
g=eye(2); 
r=l/10000; 
h=[ 1 0] ; 
cs=[ N 0]; 



X 

XXXXXXXXXXXX Determine the discrete coeficient of the system XXXXXXXXX 
X The Phi's and Gamma's of the system 
X 
% 

[ phih,gamh] =c2d(ah,bh,dt) ; 

[ phip,gamp] =c2d(ap,bp,dt) ; 

[ phir ,gamrj =c2d(ar,br ,dt) ; 

[phis,gams]=c2d(as,bs,dt) ; 

% 



% 

XXXXXXXXXXXX Generate Input for the system %%%%%%%%%%%%%%%%%%%%%%%%%% 



X 

X 

t=( 0: .1: 100)'; 
rand( 'normal' ); 
v=sqrt(r)*rand( 1000,1); 

X 

X 

us=dlsim(phis ,gams , cs ,d, v); 
X 
X 



X Time axis 

X Random are uniformly distributed 
X in the interval (0.0, 1.0). Random 
X ('normal') switches to normal with 
X mean of 0 and variance of 1 



XXXXXXXXXXXX Generate output values for the system XXXXXXJoXXXXXXXXXXXX 

X 

X 

yh=dlsim(phih,gamh,c,d,us); 
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yp=d Is im( ph ip , gamp , c , d , us ) ; 
yr=dlsim(phir ,gamr ,c,djus); 



%7 0 7o%7o%°A%7o7o%%7o Plot the generated output vs time of the system 
% 

plot(yh) ,title( ' Ships Heave Motion'); 
xlabel('Time ( sec) ') ,ylabel( ' Heave angle (rad)'); 
delete pitch, met 
meta pitch 

plot(yp) jtitle( ' Ships Pitch Motion'); 
xlabel('Time (sec) ') ,ylabel( ' Pitch angle (rad)'); 
meta 

plot(yr) ,title( ' Ships Roll Motion'); 

xlabel( 'Time (sec) ') ,ylabel( 'Roll angle (rad)'); 

meta 

plot(v) ,title( 'Random Noise vs Time'); 
xlabel('Time ( sec) ' ) ,ylabel( ' Noise Magnitude'); 
meta 

plot(us) jtitle( ' Swell with White Noise vs Time'); 
xlabel('Time (sec) ') ,ylabel( ' Swell angle (rad)'); 
meta 



%%%%%%% 
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APPENDIX B. MATLAB PROGRAM 2 

This Matlab computer program calculates the Observed Ship's Motion and esti- 
mates in real time such as the heave, pitch and roll of the ship. Required coefficients 
must be available to run the program. 



% This Matlab computer program calculates the Observed and Estimated 
% state of the Ship s Heave, Pitch and Roll based on the given data. 

% 

% 

%%%%%%%%%%%%%%%%%%%%%%%%%% SWELL parameter %%%%%%%%%%%%%%%%%%%%%%%%%%% 

% 

% 

% 

diary p. m 

% Swell period (Ts) seconds 
Ts=7 ; 



% Swell angular frequency 

ws= 2*pi*l/Ts ; 



% Wave Height (N) feet 
N=ll; 



7o Restoring coefficient 

cs=ws**2 ; 

% 

% 

7o7o7o7c7o7o7o7o7o7o%7o7o7o%7o7o7o%7o7o%7o7o7o7o Heave parameter %%%%%%%%%%%%%%%%%%%%%%%%%%% 
7o 

7o Heave period (th) seconds typical for carriers 



7o Heave angular frequency (wh) rad/sec 
wh= 2*pi*l/th ; 

7o 

7> Heave restoring cefficient (cz) 
cz=wh ; 



7> Heave damping factor (zetah) 
zetah= 0. 6; 



Damping coefficient (bh) 

bz=2*zetah*wh ; 



% 

% 

%%%%%%%%%%%%%%%%%%%%%%%%% Pitch parameter %%%%%%%%%%%%%%%%%%%%%%%%%%% 



Pitch period (tp) seconds typical for carriers 
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tp=10; 



% 

% 



Pitch angular frequency (wp) rad/sec 
wp= 2*pi*l/tp ; 



% Pitch Restoring coefficient 

cp=wp ; 

% 

% Heave damping factor (zetap) 

zetap= 0. 6 ; 

% 

% Damping coefficient (bp) 

bp=2*zetap*wp ; 

% 

% 

%%%%%%%%%%%%%%%%%%%%%%%%% Roll parameter %%%%%%%%%%%%%%%%%%%%%%%%%%X,%% 

% 

% 

% Roll period (tr) seconds typical for carriers 
tr= 18; 

% Roll angular frequency (wr) rad/sec 
wr= 2*pi*l/tr ; 



% Roll Restoring coefficient 

cr=wr ; 

% 

7 0 Roll damping factor (zetar) 
zetar= 0. 6 ; 

% Damping coefficient (br) 

br=2*zetar*wr ; 

% 

% 

%%%%%%%%%%%%%% %%%%%%%%%%% Parameters of the system %%%%%%%%%%%%%%%%%%%% 

% 

7c where: 

7c a=[];-- the A martix coefficient in continous time 

To b=[ j ; — the B matrix coefficient in continuos time 

% 

7c Swell Parameters 



% 

7c 

7c 

7c 

% 



% 



as=-0 1; -cs 0-; 
fs=[ 0; 1] ; 



Heave Parameters 



ah=[ 0 1; -cz**2 -bz] ; 
fh=[ 0; . 3] ; 
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Pitch Parameters 



diary 

% 

% 

7o 

% 

% 

% 



ap=[ 0 1; -cp**2 -bp] 
fp=[ 0; . 3] ; 



Roll Parameters 



ar=[ 0 1; -cr**2 -br] ; 
f r=[ 0; . 6] ; 
off 



Measurement Parameters 

c=[ 1 0] ; 

dt=0. 1; 

Swell Noise 
q=l/100; 

Measurement Noise 
r=l/ 8500; 



d=[ 0]; 



% 

%%%%%%%%%%%%%%%%% Determine the Discrete coeficient of the system %%%% 

% 

% Swell 

% 



% 

[phis, gams] =c2d(as,fs,dt) 

% 

% HEAVE 

% 

% 

[ phih , gamh] =c2d( ah , f h , dt ) 

% 



% 

7o 

7o 

7o 

% 



PITCH 

[ phip , gamp] =c2d( ap , f p , dt ) 
ROLL 



[phir,gamr] =c2d( ar , fr ,dt) 

% 
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% 

%%%%%%%%%% Calculate the Kalman Gains for Heave, Pitch and Roll %%%%%%% 

% 

kh=dlqe(phih,gamh,c,q,r) ; 
kp=dlqe(phip,gamp,c,q,r) ; 
kr=dlqe(phir,gamr,c,q,r) ; 

% 

%%%%%%%%%%%%%%%%%%%%%% Initial Condit ion %%%%%% %%%%%%%%%%%%%%%%% 

% 

xh=[ 0. 05; 0. 05] ; 
xhe=[ 0.05; 0.05] ; 

xp=[ 0. 05; 0. 05] ; 
xpe=[0. 05; 0. 05] ; 

xr=[ 0.050; 0.0506] ; 
xre=[ 0. 05; 0. 05] ; 
xs=[ 0. 05; 0. 05] ; 

% 

%%%%%%%%%%%%%%%%%%%%% Random Noise Generators %%%%%%%%%%%%%%%%%% 

% 

rand( 'normal' ) 
w=sqrt(q)*rand( 1 , 1500); 
v=sqrt( r)*rand( 1, 1500); 
t=0: . 1: 100; 

% 

%%%%%%%%%%%%%%%%%%%% S imu lat ion Loop %%%%%%%%%%%%%%%%%%%%%%%%% 

% 

for k=l: 1000; 

% 

%%%%%%%%%%%%%%%%%%%% Input for Heave, Pitch nad Roll %%%%%%%%% 

% 

xs(: ,k+l)=phis*xs( : ,k) + gams*w(k); 

% 

%%%%%%%%% State calculation for Heave, Pitch and Roll %%%%%%%%% 

% 

xh(: ,k+l)=phih*xh( : ,k) + gamh*xs(k) ; 
xp(: ,k+l )=phip*xp( : ,k) + gamp*xs(k) ; 
xr(: ,k+l)=phir*xr( : ,k) + gamr*xs(k) ; 

% 

% 

%%%%%%%%% Output calculation for Heave, Pitch and Roll in Discrete form 

% 

yh(: ,k+l)=c*xh(: ,k+l) + v(k+l); 
yp(: ,k+l)=c*xp(: ,k+l) + v(k+l); 
yr( : ,k+l)=c*xr(: ,k+l) + v(k+l); 

% 

% 

%%%%%%%%% State estimate calculation for Heave, Pitch and Roll %%%%%%%%%% 

% 

xhe(: ,k+l)= phih*xhe( : ,k) + kh*(yh(: ,k+l) -c*(phih*. . . 
xhe(: ,k) )); 

xpe(: ,k+l)= phip*xpe(: ,k) + kp*(yp( : ,k+l)-c*(phip*. . . 
xpe(: ,k) )); 

xre( : ,k+l)= phir*xre(: ,k) + kr*(yr(: ,k+l)-c*(phir*. . . 
xre(: ,k) )); 
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end 



% 

% 

plot(t,xh(l,: ) ,t,xhe(l,: )) ,title( 'Ships 
Observed Heave and Heave Estimate') 

xlabel( 'Time(sec) ') ,ylabel( 'Heave Angle (rad)') 
delete esth. met 
meta est 

plot(t ,xp( 1 ,: ) ,t ,xpe( 1 ,: )) ,title( ' Ships 
Observed Pitch and Pitch Estimate') 

xlabel( 'Time(sec) ') ,ylabel( ' Pitch Angle (rad)') 
meta 

plot(t ,xr( 1 , : ) ,t ,xre( 1 , : ) ) ,title( ' Ships 
Observed Roll and Roll Estimate') 

xlabel( 'Time(sec) ' ) ,ylabel( 'Roll Angle (rad)') 
meta 
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APPENDIX C. MATLAB PROGRAM 3 



This Matlab computer program calculates the Observed Ship's Motion and predicts 
in real time such as the heave, pitch and roll of the ship. Required coefficients must be 
available to run the program. 



X This Matlab computer program calculates the Observed and Predicted 
X state of the Ship s Heave, Pitch and Roll based from the given data. 

% 

X 

X 

XXXXXXXXXXXXXXXXXXXXXXXXXX SWELL parameter XXXXXXXXXXXXXXXXXXXXXXXXXXX 
X 
X 
X 

X Swell period (Ts) seconds 

Ts=7 ; 

X Swell angular frequency 

ws= 2*pi*l/Ts ; 

X Wave Height (N) feet 
N=ll| 



X 



% 



Restoring coefficient 
cs=ws**2 ; 



XXXXXXXXXXXXXXXXXXXXXXXXXX Heave parameter XXXXXXXXXXXXXXXXXXXXXXXXXXX 
X 

X Heave period (th) seconds typical for carriers 
th= 11 ; 

X 

X Heave angular frequency (wh) rad/sec 
wh= 2*pi*l/th ; 

X 

X Heave restoring cefficient (c z) 

cz=wh ; 

X Heave damping factor (zetah) 

zetah= 0. 6; 

% Damping coefficient (bh) 

bz=2*zetah*wh ; 

% 

% 

XXXXXXXXXXXXXXXXXXXXXXXXX Pitch parameter XXXXXXXXXXXXXXXXXXXXXXXXXXX 

X 

X 

X Pitch period (tp) seconds typical for carriers 
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tp=10; 



% Pitch angular frequency (wp) rad/sec 

wp= 2*pi*l/tp ; 

% 

% 

% Pitch Restoring coefficient 



% Heave damping factor (zetap) 

zetap= 0. 6 ; 

% 

% Damping coefficient (bp) 

bp=2*zetap*wp ; 

% 

% 

%%%%%%%%%% %%%%%%%%%%%%%%% Roll parameter %%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% 

% 

% Roll period (tr) seconds typical for carriers 
tr= 18; 

% Roll angular frequency (wr) rad/sec 
wr= 2*pi*l/tr ; 



% Roll Restoring coefficient 

cr=wr ; 

% 

7o Roll damping factor (zetar) 

zetar= 0. 6 ; 

% Damping coefficient (br) 

br=2*zetar ,v wr ; 

% 

7o 

%%%%%%%%%%%%%% %%%%%%%%%%% Parameters of the system %%%%%%%%%%%%%%%%%%%% 

% 

% where: 

% a=[ ] ; — the A martix coefficient in continous time 

% b=[ ] ; — the B matrix coefficient in continuos time 

% 

% Swell Parameters 

% ======= = ====^== 

% 

as=[ 0 1; -cs 0] ; 
fs=[ 0; 1] ; 

7a 

% 

% Heave Parameters 

% 

ah=[ 0 1; -cz**2 -bz] ; 
fh=[ 0; . 3] ; 

% 
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Pitch Parameters 



% 

% 



% 

% 

7o 

7o 

% 



% 

% 

7o 

7o 

7o 

% 



7o 

7, 

7o 

7 0 

7o 

7o 



ap=[ 0 1; -cp**2 -bp] ; 
fp=[ 0; . 3] j 



Roll Parameters 



ar=[ 0 l;-cr**2 -br] ; 
f r=[ 0; . 6] ; 



Measurement Parameters 



c=[ 1 0] ; 
dt=0. 1; 

Swell Noise 

q=l/ 100; 

Measurement Noise 
r=l/3500; 



d=[ 0]; 

% 

% 

7o7o7c7o7o7 a 7o7o7 a 7o7o7o7 a 7o7o7o7o Determine the Discrete coeficient of the system %%%% 



7o 

7o Swell 

7o 



7o 

7o 

% 

% 

% 

7o 

7o 

7o 

7o 

7o 

7o 

% 

7o 



[ phis , gams] =c2d( as , fs ,dt) 
HEAVE 

[ phih,gamh] =c2d(ah, fh,dt) 
PITCH 

[ phip , gamp] =c2d( ap , f p , dt ) 
ROLL 

[ phir ,gamr] =c2d(ar, fr,dt) ; 
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% 

%%%%%%%%%% Calculate the Kalman Gains for Heave, Pitch and Roll %%%%%%% 

% 

kh=dlqe(phih,gamh, c,q,r) ; 
kp=dlqe(phip,gamp,c,q,r) ; 
kr=dlqe(phir ,gamr ,c,q, r) ; 

% 

%%%%%%%%%%%%%%%%%%%%%% Initial Condition %%%%%%%%%%%%%%%%%%%%%%% 

% 

xh=[ -0. 05; -0. 05] ; 
xhe=[ -0. 05; -0. 05] ; 

xp=[0. 05; -0. 05] ; 
xpe=[ -0. 05; -0. 05] ; 

xr=[0. 050; -0. 050] ; 
xre=[ -0. 05; -0. 05] ; 

xs=[0. 05; -0. 05] ; ; 

% 

%%%%% %%%%%%%%%%%%%%%% Random Noise Generators %%%%%%%%%%%%%%%%%% 

% 

rand( 'normal* ) 
w=sqrt(q)*rand( 1,1200); 
v=sqrt( r)*rand( 1 , 1200); 
t=0: . 1: 95; 

% 

%%%%%%%%%%% %%%%%%%%% Simulation Loop %%%%%%%%%%%%%%%%%%%%%%%%% 

% 

for k=l: 950; 

% 

%%%%%%%% %%%%%%%%%%%% Input for Heave, Pitch nad Roll %%%%%%%%% 

% 

xs(: ,k+l)=phis*xs( : ,k) + gams*w(k); 

%%%%%%%%% State calculation for Heave, Pitch and Roll %%%%%%%%% 

% 

xh(: ,k+l)=phih*xh( : ,k) + gamh*xs(k) ; 
xp(: ,k+l)=phip*xp( : ,k) + gamp*xs(k) ; 
xr(: ,k+l)=phir*xr( : ,k) + gamr*xs(k) ; 

% 

% 

%%%%%% Output calculation for Heave, Pitch and Roll in Discrete form 

yh(: ,k+l)=c*xh(: ,k+l) + v(k+l); 
yp(: ,k+l)=c*xp(: ,k+l) + v(k+l); 
yr(: ,k+l)=c*xr(: ,k+l) + v(k+l); 

% 

% 

7o7o7o°/o7o7o7o State estimate calculation for Heave, Pitch and Roll %%%%%%%%%% 

% 

xhe(: ,k+l)= phih*xhe( : ,k) + kh*(yh( : ,k+l) -c*(phih*. . . 
xhe(: ,k) )); 

xpe(: ,k+l)= phip*xpe(: ,k) + kp*(yp(: ,k+l)-c*(phip*. . . 
xpe(: ,k) )); 

xre(:,k+l)= phir*xre( : ,k) + kr*(yr( : ,k+l) -c*(phir*. . . 
xre(: ,k) )); 



% 

%7o%%%% State Prediction Calculation for Heave, Pitch and Roll %%%%%%% 
% 

xhp(: ,k+l)=phih*xhe( : ,k) + gamh*xs(k) ; 
xpp(: ,k+l)=phip*xpe(: ,k) + gamp*xs(k) ; 
xrp(: ,k+l)=phir*xre( : ,k) + gamr ,v xs(k) ; 

end 

% 

% 



plot(t,xh(l,: ,t,xhp(l,: )),title( 'Ships 

Observed Heave and Predicted Heave') 

xlabel( 'Time(sec) ') ,ylabel( 'Heave Angle (rad)') 
delete predicth. met 
meta predicth 

plot(t ,xp( 1 ,:),'- ' , t ,xpp( 1, : )) , title( ' Ships 
Observed Pitch and Predicted Pitch' ) 

xlabel( 'Time(sec)' ) ,ylabel( 'Pitch Angle (rad)') 
meta 

plot(t ,xr( 1 , : ) , ' -' ,t ,xrp( 1, : ) ) ,title( ' Ships 
Observed Roll and Predicted Roll') 

xlabel( 'Time(sec) ' ) ,ylabel( 'Roll Angle (rad)') 
meta 
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